{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "839a3961",
   "metadata": {},
   "source": [
    "## Example: 3D Heat Equation\n",
    "\n",
    "### What Features Are We Using\n",
    "\n",
    "* Mesh data \n",
    "* Domain Decomposition\n",
    "\n",
    "We now look at a more complicated example at and show how simulation results\n",
    "can be visualized. This example solves the heat equation,\n",
    "\n",
    "$$\\frac{\\partial\\phi}{\\partial t} = \\nabla^2\\phi$$\n",
    "\n",
    "using forward Euler temporal integration on a periodic domain.  We could use a\n",
    "5-point (in 2D) or 7-point (in 3D) stencil, but for demonstration purposes we\n",
    "spatially discretize the PDE by first constructing (negative) fluxes on cell faces, e.g.,\n",
    "\n",
    "$$F_{i+^1\\!/_2,\\,j} = \\frac{\\phi_{i+1,j}-\\phi_{i,j}}{\\Delta x}$$\n",
    "\n",
    "and then taking the divergence to update the cells:\n",
    "\n",
    "   $$\\phi_{i,\\,j}^{n+1} = \\phi_{i,\\,j}^n\n",
    "   + \\frac{\\Delta t}{\\Delta x}\\left(F_{i+^1\\!/_2,\\,j}-F_{i-^1\\!/_2,\\,j}\\right)\n",
    "   + \\frac{\\Delta t}{\\Delta y}\\left(F_{i,\\,j+^1\\!/_2}-F_{i,\\,j-^1\\!/_2}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08b4ec5f",
   "metadata": {},
   "source": [
    "The code to generate the initial condition is in `mykernel.H` and looks like: \n",
    "```C++\n",
    "{\n",
    "    Real x = prob_lo[0] + (i+Real(0.5)) * dx[0];\n",
    "    Real y = prob_lo[1] + (j+Real(0.5)) * dx[1];\n",
    "    Real z = prob_lo[2] + (k+Real(0.5)) * dx[2];\n",
    "    Real r2 = ((x-Real(0.25))*(x-Real(0.25))+(y-Real(0.25))*(y-Real(0.25))+(z-Real(0.25))*(z-Real(0.25)))/Real(0.01);\n",
    "    phi(i,j,k) = Real(1.) + std::exp(-r2);\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33bd71f2",
   "metadata": {},
   "source": [
    "## Running the code\n",
    "\n",
    "The simulation can be ran as `./03_HeatEquation inputs`. \n",
    "\n",
    "The following inputs parameters could be tweaked:\n",
    "\n",
    "```\n",
    "nsteps        = 1000              # number of time steps to take\n",
    "plot_int      = 100               # write plots every n steps\n",
    "n_cell        = 128               # number of cells in the domain\n",
    "max_grid_size = 64                # max grid size used for domain decomposition\n",
    "\n",
    "```\n",
    "\n",
    "Although we are running this example in serial, we decompose the domain into multiple boxes, anticipating more complicated problems where we have mesh refinement:\n",
    "\n",
    "```C++\n",
    "        IntVect dom_lo(AMREX_D_DECL(       0,        0,        0));\n",
    "        IntVect dom_hi(AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1));\n",
    "        Box domain(dom_lo, dom_hi);\n",
    "\n",
    "        // Initialize the boxarray \"ba\" from the single box \"bx\"\n",
    "        ba.define(domain);\n",
    "        // Break up boxarray \"ba\" into chunks no larger than \"max_grid_size\" along a direction\n",
    "        ba.maxSize(max_grid_size);\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b772b897",
   "metadata": {},
   "source": [
    "## Visualizating the results\n",
    "\n",
    "Below we give some python code to visualizate the solution using yt:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29609c78",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yt\n",
    "from yt.frontends import boxlib\n",
    "from yt.frontends.boxlib.data_structures import AMReXDataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18ec883c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb533aa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = AMReXDataset(\"plt00000\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d834eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.field_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "899767fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "sl = yt.SlicePlot(ds, 2, ('boxlib', 'phi'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "679eaca7",
   "metadata": {},
   "outputs": [],
   "source": [
    "sl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72af42b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = AMReXDataset(\"plt01000\")\n",
    "sl = yt.SlicePlot(ds, 2, ('boxlib', 'phi'))\n",
    "sl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd7332dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "sl.annotate_grids()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
